 subroutine getpqr_ri(icore, s_in, irreps_in, ri1,ri2, popp, popq, ioppq, &
                      popr, pops, ioprs, ioppr, isym)
! One should note that the memory in getpqr is nabc*2 
! so the memory size of icore must be carefully treatted
! ispin = 1 aa; ispin = 2 bb; ispin=3 ab
 use mod_ri
 use mod_iop
 use mod_symm
 use mod_orbit
 implicit none 
!----------------------------------
 integer, intent(in) ::  s_in, irreps_in,ioppq, ioprs, isym, ioppr
 real*4, intent(in) :: ri1(*), ri2(*)
 integer, intent(in) :: popp(8),popq(8),popr(8)
 integer, intent(in) :: pops(8)
 real*4, intent(out) :: icore(*)

 integer irrep, n, npqr, npqr2
 integer ioff,  ioffout
 integer i0, i1, s, irreps
 integer nsize,  npqr2
 integer ioffvo(8)
 integer, external :: irpdso, idsymsoc, isymoffso
 real*8, external :: nrm2so
   
    zorb(1:8) = 0 
    zorb(1) = 1
          irreps = irreps_in
          s = s_in
          npqr = idsymsoc(irreps, popp, popq, 1, popr, zorb, 1)
          i0 = 1
          i1 = i0 + max(npqr*isd, idsymsoc(irreps, naux8, zorb, 1, popr, zorb, 1))
          call getpqrso_incore(icore, s, irreps, naux8, zorb, 1, popr, pops, ioprs, ri2, 1)
          call VmnpqVmnrs_to_Vpqrs(ri1, icore, 1, irreps, naux8, zorb, 1, popp, popq, ioppq, &
                        popr, zorb, 1, icore(i1), irreps, 1.d0, 0.d0)
          if(ioppq==2) then
             call iexpso(icore(i1), icore(i0), popp, popq, 2, 1, popr, zorb, 1, 0, irreps)
             npqr2 = idsymsoc(irreps, popp, popq, 1, popr, zorb, 1)
             call myicopyso(icore(i0), icore(i1), npqr2)
          endif
          !(pq|rs) to <pr|qs> with prq
          call sstgenso(icore(i1), icore,  nsize, popp, popq, popr, zorb, icore,&
                        irreps, '1324')
!!!!<pr||qs> = <pr|qs> - <rp|qs>
          if(ioppr==0) then 
             npqr = idsymsoc(irreps, popp, popr, 1, popq, zorb, 1)
             call isymtrso(icore, icore(i1), popp, popr, 1, 1, popq, zorb, 1, 0, irreps)
             call axpyso(npqr, -1.d0, icore(i1), 1, icore, 1)
             call icompso(icore, icore(i1), popp, popr, 1, 0, popq, zorb, 1, 1,irreps) 
             npqr = idsymsoc(irreps, popp, popr, 0, popq, zorb, 1) 
             call myicopyso(icore(i1), icore, npqr)
          elseif(ioppr==2) then 
             npqr = idsymsoc(irreps, popp, popr, 1, popq, zorb, 1)
             call isymtrso(icore, icore(i1), popp, popr, 1, 1, popq, zorb, 1, 0, irreps)
             call axpyso(npqr, 1.d0, icore(i1), 1, icore, 1)
             call icompso(icore, icore(i1), popp, popr, 1, 2, popq, zorb, 1, 1,irreps) 
             npqr = idsymsoc(irreps, popp, popr, 2, popq, zorb, 1) 
             call myicopyso(icore(i1), icore, npqr)
          endif
  return
  end
